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Abstract. The fractal or Hausdorff dimension is a measure of rough- 
ness (or smoothness) for time series and spatial data. The graph of 
a smooth, differentiable surface indexed in M*^ has topological and frac- 
tal dimension d. If the surface is nondifferentiable and rough, the fractal 
dimension takes values between the topological dimension, d, and d+1. 
We review and assess estimators of fractal dimension by their large sam- 
ple behavior under infill asymptotics, in extensive finite sample simula- 
tion studies, and in a data example on arctic sea-ice profiles. For time 
series or line transect data, box-count, Hall-Wood, semi-periodogram, 
discrete cosine transform and wavelet estimators are studied along with 
variation estimators with power indices 2 (variogram) and 1 (mado- 
gram), all implemented in the R package f ractaldim. Considering both 
efficiency and robustness, we recommend the use of the madogram es- 
timator, which can be interpreted as a statistically more efficient ver- 
sion of the Hall- Wood estimator. For two-dimensional lattice data, we 
propose robust transect estimators that use the median of variation 
estimates along rows and columns. Generally, the link between power 
variations of index p > for stochastic processes, and the Hausdorff 
dimension of their sample paths, appears to be particularly robust and 
inclusive when p=l. 

Key words and phrases: Box-count, Gaussian process, Hausdorff di- 
mension, madogram, power variation, robustness, sea-ice thickness, 
smoothness, variogram, variation estimator. 
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1. INTRODUCTION 

Fractal-based analyses of time series, transects, 
and natural or man-made surfaces have found exten- 
sive applications in almost all scientific disciplines 
(Mandelbrot, 1982). While much of the literature 
ties fractal properties to statistical self-similarity, no 
such link is necessary. Rather, we adopt the argu- 
ment of Bruno and Raspa (1989), Davies and Hall 
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(1999) and Gneiting and Schlather (2004) that the 
fractal or Hausdorff dimension quantifies the rough- 
ness or smoothness of time series and spatial data 
in the limit as the observational scale becomes in- 
finitesimally fine. In practice, measurements can only 
be taken at a finite range of scales, and usable esti- 
mates of the fractal dimension depend on the avail- 
ability of observations at sufficiently fine temporal or 
spatial resolution (Malcai et al., 1997; Halley et al., 
2004). 

We follow common practice in defining the fractal 
dimension of a point set X C M'^ to be the classi- 
cal Hausdorff dimension (Hausdorff, 1919; Falconer, 
1990). For e > 0, an e-cover of X is a finite or count- 
able collection {Bi -.i = 1,2, .. .} of balls Bi C M'^ of 
diameter \Bi\ less than or equal to e that covers X. 
With 



H%X) 



lim inf < 

£^0 



f2\Bif-{Bi:i = 1,2,...} 



is an e-cover of X j 

denoting the ^-dimensional Hausdorff measure of X , 
there exists a unique nonnegative value D such that 
H^X) = oo if (5 < and H^{X) = 0if6>D. This 
value D is the Hausdorff dimension of the point 
set X. Under weak regularity conditions, the Haus- 
dorff dimension coincides with the box-count dimen- 
sion, 

logiV(e) 



1 



Dbc 



■ lim , / , N 
e-^o log(l/e) 



where N{e) denotes the smallest number of cubes of 
width e in M.'^ which can cover X, and also with other 
natural and/or time- honored notions of dimension 
(Falconer, 1990). 

In this paper we restrict attention to the case in 
which the point set 



x = {{t,Xt)e 



is the graph of time series or spatial data observed 
at a finite set T C M'' of typically regularly spaced 
times or locations. The fractal dimension then refers 
to the properties of the curve {d = 1) or surface 
(d > 2) that arises in the continuum limit as the 
data are observed at an infinitesimally dense subset 
of the temporal or spatial domain, which without 
loss of generality can be assumed to be the unit in- 
terval or unit cube. In time series analysis and spa- 
tial statistics, this limiting scenario is referred to as 
infill asymptotics (Hall and Wood, 1993; Dahlhaus, 
1997; Stein, 1999). 



If the limit curve or limit surface is smooth and 
differentiable, its fractal dimension, D, equals its 
topological dimension, d. For a rough and nondif- 
ferentiable curve or surface, the fractal dimension 
may exceed the topological dimension. For example, 
suppose that {Xt : t G W^} is a Gaussian process with 
stationary increments, whose variogram or structure 
function, 

(2) j2{t) = ^E{Xu-Xu+tf, 
satisfies 

(3) 72(t) = |c2tr + 0(|ir+^) ast^O, 

where a G (0, 2], /3 > 0, and C2 > 0, and | • | denotes 
the Euclidean norm. Then the graph of a sample 
path has fractal dimension 

(4) D = d+l-^ 

almost surely (Orey, 1970; Adler 1981). This rela- 
tionship links the fractal dimension of the sample 
paths to the behavior of the variogram or structure 
function at the coordinate origin, and can be ex- 
tended to broad classes of potentially anisotropic 
and nonstationary processes, as well as some non- 
Gaussian processes (Adler, 1981; Hall and Roy, 1994; 
Xue and Xiao, 2011). It allows us to think of fractal 
dimension as a second-order property of a Gaussian 
stochastic process, in addition to being a roughness 
measure for a realized curve or surface. Accordingly, 
we refer to the index a in the asymptotic relation- 
ship (3) as the fractal index. 

Table 1 provides examples of Gaussian processes 
that exhibit this asymptotic behavior. Fractional 
Brownian motion is a nonstationary, statistically self- 
similar process that is defined in terms of the vari- 
ogram (Mandelbrot and Van Ness, 1968). The other 
entries in the table refer to stationary processes with 
covariance function a{t) = cov{Xu,Xt-\.u), which re- 
lates to the variogram as 

72(t) =cj(0) -cj(t), t£R'^. 

Key examples include the Matern family (Matern, 
1986; Guttorp and Gneiting, 2006), used by Goff 
and Jordan (1988) to parameterize the fractal di- 
mension of oceanic features; the Cauchy class, intro- 
duced by Gneiting and Schlather (2004) to illustrate 
local and global properties of random functions; and 
the Dagum family (Berg, Mateu and Porcu, 2008). 
In simulation settings, the powered exponential class 
(Yaglom, 1987) is a convenient default example. The 
exponent /3 in the asymptotic relationship (3) equals 
/3 = 2 — a for the Matern class, /3 = a for the powered 
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Table 1 

Some parametric classes of variograms and covariance functions for a Gaussian 
process {Xt :t G R''}. The covariance functions have been normalized such that 
(7(0) — 1. Here, a is the fractal index, c> is a range parameter, and K^, is a modified 
Bessel function of the second kind of order v . The Matern and Dagum families allow 
for less restrictive assumptions on the parameters than stated here 

Class Variogram or covariance Parameters 



Fractional Brownian motion 72 (t) = \ct\°' a £ (0, 2] 

Matern G{t) = ^^(1/';/ M"'^ K^/^{\ct\) a e (0, 2] 

Powered exponential a{t) = &cp[—\ct\°') afE(0,2] 

Cauchy cr(f) = (l + |ciP)-"/" aG(0,2];r>0 

Dagum <^(t) = i_(_Ml_)-A r G (0, 2]; a e (0, r) 



exponential and Cauchy families, and (3 = t for the 
Dagum family. Generally, the smaller the value of /3, 
the harder the estimation of the fractal index, a, 
and the more pronounced the finite sample bias of 
estimators of the fractal index or fractal dimension. 

As an illustration for time series or line transect 
data. Figure 1 displays Gaussian sample paths from 
the powered exponential family with the fractal in- 
dex, a, ranging from 1.9 to 0.2, and the fractal di- 
mension, D, extending from 1.05 to 1.9. The small- 
er a, the rougher the sample path, and the larger the 
fractal dimension, with the lower limit, D = 1, being 
associated with a smooth curve, and the upper limit, 
D = 2, corresponding to a space-filling, exceedingly 
rough graph. Figure 2 shows a realization from the 
nonstationary Gaussian Matern model of Anderes 
and Stein (2011), in a case in which the fractal in- 
dex and the fractal dimension vary linearly along the 
unit interval. To illustrate the visual effects of the 
measurement scale, both the original sample path of 
size 10,000 and an equidistantly thinned version of 
size 1,000 are shown. 

Turning to spatial data, Figure 3 shows Gaussian 
sample surfaces from the powered exponential class 
with the fractal index, a, being equal to 1.5, 1.0 
and 0.2. The surfaces are increasingly rough with 
the fractal dimension, D, being equal to 2.25, 2.5 
and 2.9, respectively. 

A wealth of applications requires the characteri- 
zation of the roughness or smoothness of time series, 
line transect or spatial data, with Burrough (1981) 
and Malcai et al. (1997) summarizing an impressive 
range of experimental results. For example, fractal 
dimensions have been studied for geographic pro- 
files and surfaces, such as the underside of sea ice 
(Rothrock and Thorndike, 1980), the topography 
of the sea floor (Goff and Jordan, 1988), Martial 



surface (Orosei et al., 2003) and terrestrial features 
(Weissel, Pratson and Malinverno, 1994; Turcotte, 
1997; Gagnon, Lovejoy and Schertzer, 2006). Fur- 
ther references can be found in Molz, Rajaram and 
Lu (2004) for applications in subsurface hydrology 
and in Halley et al. (2004) for applications in ecol- 
ogy. Not surprisingly then, estimators for the frac- 
tal dimension have been proposed and widely used 
in various literatures, including physics, engineer- 
ing, the earth sciences, and statistics, with the works 
of Burrough (1981), Goff and Jordan (1988), Bruno 
and Raspa (1989), Dubuc et al. (1989a, 1989b), Jake- 
man and Jordan (1990), Theiler (1990), Klinkenberg 
and Goodchild (1992), Schepers, van Beek and Bass- 
ingthwaighte (1992), HaU and Wood (1993), Con- 
stantine and Hall (1994), Kent and Wood (1997), 
Davies and Hah (1999), Chan and Wood (2000), Zhu 
and Stein (2002), Chan and Wood (2004) and Bez 
and Bertrand (2011) being examples. Our objectives 
in this paper are to survey the literature across disci- 
plines, review and assess the various types of estima- 
tors, and provide recommendations for practition- 
ers, along with novel directions for theoretical work. 

The remainder of the paper is organized as follows. 
In Section 2 we describe estimators of fractal dimen- 
sion for time series and line transect data, including 
box-count, Hall-Wood, variogram, madogram, pow- 
er variation, semi-periodogram and wavelet-based 
techniques. Then in Section 3 we assess and compare 
the estimators. Considering both efficiency and ro- 
bustness, we concur with Bruno and Raspa (1989) 
and Bez and Bertrand (2011) and recommend the 
use of the madogram estimator, that is, the varia- 
tion estimator with power index p = l, which can be 
interpreted statistically efficient version of the 
Hall- Wood estimator. An underlying motivation is 
that for an intrinsically stationary process with var- 
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Fig. 1. Sample paths of stationary Gaussian processes with powered exponential covanance function, = exp( — and 
fractal index, ot, equal to 1.9, I.4, 1-0 and 0.2. The corresponding values of the fractal dimension, D, are 1.05, 1.3, 1.5 and 
1.9, respectively. The simulation domain is a grid over the unit interval with spacing 1/1,024. 



I. 
I 




1 



4 - 



*t - 




— r — 



— 1 — 
01 



— t— 



— 1 — 



04 



Fig. 2. A sample path from the nonstationary Gaussian Matern process of Anderes and Stem (2011), where the fractal 
dimension, D, grows linearly from D = 1 at time t = to D — 2 at time t — 1. To illustrate the visual effects of the temporal 
resolution, both the original sample path with grid spacing 1/10,000 (top panel) and an equidistantly thinned version with 
grid spacing 1/1,000 (bottom panel) are shown. The nonstationary covariance is given by equation (10) of Anderes and Stein 
(2011) with — 1, p= 1/2, and linearly varying local smoothness parameter, ut = 1 ~t. 
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Fig. 3. Sample surfaces of stationary spatial Gaussian processes with the powered exponential covariance function, 
(j(t) — exp{ — \t\") , and fractal index, a, equal to 1.5, 1.0 and 0.2. The corresponding values of the fractal dimension, D, 
are 2.25, 2.5 and 2.9. The simulation domain is a grid on the unit square m with spacing 1/512 along each coordinate. 



iogram of order p > of the form 



(5) 



■X 



u+t\ 



the relationship (4) between the fractal index, a, 
and the fractal dimension, D, appears to be more 
robust and inclusive when p = 1, as compared to 
the default case, in which p = 2. 

Section 4 discusses ways in which estimators for 
the time series or line transect case can be adapted 
to spatial data observed over a regular lattice in M? , 
and Section 5 evaluates these proposals. Our pre- 
ferred estimators in this setting are simple, robust 
transect estimators that employ the median of varia- 
tion estimates with power index p = l along individ- 
ual rows and columns. A data example on arctic sea- 
ice profiles is given in Section 6. The paper ends in 
Section 7, where we make a call for new directions in 
theoretical and methodological work that addresses 
both probabilists and statisticians. Furthermore, we 
hint at inference for nonstationary or multifractional 
processes, where the fractal dimension of a sample 
path may vary temporally or spatially. All compu- 
tations in the paper use the fractaldim package 
(Sevcikova, Gneiting and Percival, 2011), which im- 
plements our proposals in R (Ihaka and Gentleman, 
1996). 

2. ESTIMATING THE FRACTAL DIMENSION 
OF TIME SERIES AND LINE TRANSECT 
DATA 

Spurred and inspired by the now classical essay 
of Mandelbrot (1982), a large number of methods 
have been developed for estimating fractal dimen- 
sion. By the early 1990s a sizable, mostly heuris- 



tic literature on the estimation of fractal dimension 
for time series and line transect data had accumu- 
lated in the physical, engineering and earth sciences, 
where various reviews are available (Dubuc et al., 
1989a; Klinkenberg and Goodchild, 1992; Schepers, 
van Beek and Bassingthwaighte, 1992; Gallant et al., 
1994; Khnkenberg, 1994; Schmittbuhl, Vilotte and 
Roux, 1995). These developments prompted the sta- 
tistical community to introduce new methodology, 
along with asymptotic theory for box-count (Hall 
and Wood, 1993), variogram (Constantine and Hall, 
1994; Kent and Wood, 1997), level crossing (Feuer- 
verger. Hall and Wood, 1994) and spectral (Chan, 
Hall and Poskitt, 1995) estimators, among others. 

Essentially all methods follow a common scheme, 
in that: 

(a) a certain numerical property, say Q, of the 
time series or line transect data is computed as a func- 
tion of "scale," say e; 

(b) an asymptotic power law Q{e) oc as the 
scale e — >■ becomes infinitesimally small is derived 
or postulated; where 

(c) the scaling exponent, 5, is a linear function of 
the fractal dimension, D; 

(d) and thus D is estimated by linear regression 
of log (5(e) on loge, with emphasis on the smallest 
observed values of the scale e. 

Table 2 shows the data property, the measure of 
scale and the scaling law for various methods. For 
techniques working in the spectral domain, the scal- 
ing law applies as the frequency grows to infinity, 
equivalent to the scale becoming infinitesimally small 
in the time domain. 

In the balance of this section, we describe the 
most popular estimators of fractal dimension in the 
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Table 2 

Some extant methods for estimating the fractal dimension of time series and line transect data 



Method 


Property 


Scale 


Scaling law 


Regime 


Box-count 


N{e): number of boxes 


e: box width 


iV(e) oc£-" 


e^O 


Divider 


Lie): length of curve 


e: step size 




e^O 


Level crossing 


M{h): number of crossings 


h: bandwidth 


M{h)(xh^-'^ 


h^Q 


Variogram 


72 (t): variogram 


t: lag 


72Woct^-^° 


t^Q 


Madogram 


71 (t): madogram 


t: lag 




t^Q 


Spectral 


/(oj): spectral density 


uj: frequency 




— j> 00 


Wavelet 


v'^ij): wavelet variance 


r: scale 




r-S>0 



equally spaced time series or line transect setting. 
Without loss of generality, we may assume that the 
observation domain is the unit interval. In the case 
of = n + 1 equally spaced observations, the data 
graph is the point set 

(6) |(t,Xt):t = ^,i = 0,l,...,n| CM^. 

The relevant asymptotic regime then is infill asymp- 
totics, in which the number of observations grows to 
infinity, whereas the underlying domain, namely the 
unit interval, remains fixed. For convenience in what 
follows, we refer to both n and Ug as the sample size. 

2.1 Box-Count Estimator 

The popular box-count estimator is motivated by 
the scaling law (1) that defines the box-count di- 
mension. The basic idea is simple, in that the time 
series graph is initially covered by a single box. The 
box is divided into four quadrants, and the num- 
ber of cells required to cover the curve is counted. 
Then each subsequent quadrant is divided into four 
subquadrants, and one continues doing so until the 
box width equals the resolution of the data, keep- 
ing track of the number of quadrants required to 
cover the graph at each step. If A'^(e) denotes the 
number of boxes required at width or scale e, the 
box-count estimator equals the slope in an ordinary 
least squares regression fit of logiV(e) on loge. Sim- 
ilarly to Mandelbrot's (1967) divider technique, the 
method can be used to quantify the fractal dimen- 
sion of any planar point set, rather than just equally 
spaced time series or line transect data. 

In our setting of a data graph of the form (6), 
where, for simplicity, we assume that n = 2^ is a pow- 
er of 2, the box-count algorithm can be summarized 
as follows. Let n = maxo<j<n X^y^ — mino<j<nXj/„ 
denote the range of the data. Consider scales = 
2^^"^ where A; = 0, 1, . . . , i^T. At the largest scale, 



e/^ = 1, the graph (6) can be covered by a single 
box of width 1 and height u, which we now call the 
bounding box. At scale the bounding box can 
be tiled by 4-^"^^ boxes of width 2''~^ and height 
u2^~^ each, and we define N{£]^) to be the number 
of such boxes that intersect the linearly interpolated 
data graph. Figure 4 provides an illustration on two 
of the datasets in Figure 1, for which n = 1024 and 
K = 10. For example, the upper left plot consid- 
ers = 8, where = 2~^ and A^(e8) = H, and the 
middle left plot looks at k = 5, where £5 = 2~^ and 
iV(e5) = 207. The naive box-count estimator then 
uses the slope in an ordinary least squares regres- 
sion fit of logA^(e) on loge, that is, 

Dec = -iY.^Sk-s)logN{ek)\iY.^Sk-sA , 

lfc=0 J U=o J 

where Sk = logefc and s is the mean of sq, si, . . . , sk- 
In our illustrating example, this leads to the esti- 
mates shown in the lower row of Figure 4. 

Several authors identified problems with the naive 
estimator that includes all scales in the regression fit 
of log A^(e) on loge, and proposed modifications that 
address these issues (Dubuc et al., 1989a; Liebovitch 
and Toth, 1989; Block, von Bloh and Schellnhuber, 
1990; Taylor and Taylor, 1991). Indeed, one always 
has A^(eo) ^n- and A^(e/<) = 1, which suggests that 
the counts at both the smallest and the largest scales 
ought to be discarded. Liebovitch and Toth (1989) 
proposed to exclude the smallest scales Sk for which 
N{£k) > n/5, as well as the two largest scales, from 
the regression fit. We adopt this proposal in our 
standard version of the box-count estimator, as il- 
lustrated in Figure 5. The restriction on the scales 
improves the statistical and computational efficiency 
of the estimator. However, it is in the limit as e — t- 
that the underlying scaling law (1) operates, and 
thus it is unfortunate that information at the very 
smallest scales is discarded. 
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Fig. 4. Illustration of the box-count algorithm and naive box-count estimates for the two datasets in the lower row of Figure 1. 
See the text for details. 



A natural variant of the box-count estimator uses 
scales £i = l/n, rather than scales = at the 
powers of 2 only. In the next section we discuss a re- 



lated estimator that takes up this idea, addresses the 
aforementioned limitations, and is tailored to time 
series data of the form (6). 



«E»1 (□■•t.Bi; a-«.I (D>1JB]| 




T 1 r 1 1 1 r r \ 1 1 1 1 1 1 r 



Fig. 5. Log-log regression for the standard version of the box-count estimator and the datasets in the lower row of Figure 1. 
Only the points marked with filled circles are used when fitting the regression line. 
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a = 1 (D = 1.5) a = 0.2 (D = 1.9) 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Fig. 6. Illustration of the Hall-Wood algorithm for the datasets in the lower row of Figure 1. The quantity A{l/n) is computed 
as the sum of the colored areas, where n — 1,024, I = 10 (upper row) and I = 30 (lower row). 



2.2 Hall-Wood Estimator 

Hall and Wood (1993) introduced a version of 
the box-count estimator that operates at the small- 
est observed scales and avoids the need for rules of 
thumb in its implementation. 

To motivate their proposal, let A{e) denote the to- 
tal area of the boxes at scale e that intersect with the 
linearly interpolated data graph (6). There are N{e) 
such boxes, and so A{e) oc N{e)e'^, which leads us to 
a reformulation of definition (1), namely, 

logA{e) 



(7) 

At scale ei 
A{l/n) is 



Dbc = 
l/n, where / 



Um 

£^0 log(e) 



1,2, 



an estimator of 



(8) 



Ail/n) 



^ Vn/i\ 



n 



i=l 



^il/n — ^(i-l)l/n\i 



where [n/l\ denotes the integer part of n/l. Fig- 
ure 6 suggests a natural interpretation of this quan- 
tity, in that it approximates A{l/n), with all features 
at scales less than l/n being ignored. For an alter- 
native, and potentially preferable, interpretation in 
terms of power variations, see Section 2.4. 



The Hall-Wood estimator with design parame- 
ter m = 1, as used in the numerical experiments of 
Hall and Wood (1993), is based on an ordinary least 
squares regression fit of log A(l/n) on log(Z/n): 



(si-s) log A{l/n 




Dnw = 2 - 
(9) 

where L>2, si = log{l/n) and s = J2i=i ^i- ^^^^ 
and Wood (1993) recommended the use of L = 2 
to minimize bias, which is unsurprising, in view of 
the limit in (7) being taken as e — )• 0. Using L = 
2 yields our standard implementation of the Hall- 
Wood estimator, namely, 



(10) Dnw = 2 



logyl(2/n) - log^(l/n) 
b^2 



Figure 7 shows the corresponding log-log plots and 
regression fits in our illustrating example. 

2.3 Variogram Estimator 

Owing to its intuitive appeal and ease of imple- 
mentation, the variogram estimator has been very 
popular. A prominent early application is that of 
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Fig. 7. Log-log regression for the Hall-Wood estimator (10) and the datasets in the lower row of Figure 1. Only the two 
points marked with filled circles are used when fitting the regression line. 



Burrough (1981). The first asymptotic study under 
the infill scenario was published in the physics liter- 
ature (Jakeman and Jordan, 1990), followed by key 
contributions of Constantine and Hall (1994), Kent 
and Wood (1997), Davies and Hah (1999) and Chan 
and Wood (2004) in statistical journals. 

Recall that the variogram or structure function 7(t) 
of a stochastic process {Xt : t G M} with stationary 
increments is defined in (2) as one-half times the ex- 
pectation of the square of an increment at lag t. The 
classical method of moments estimator for 7(t) at 
lag t = l/n from time series or line transect data (6) 
is 

1 " 

(11) y2(//n) = — — -^(X,/„-X(,_;)/„)l 

' i=l 

In view of the relationship (4) between the fractal 
index, defined in (3), and the fractal dimension, D, 
a regression fit of \ogV{t) on logt yields the vari- 
ogram estimator, 

Sv;2 = 2 - i I ^(s, - s) log % (l/n) I 



where L > 2, si = log(Z/n) and s is the mean of 
si,. . . ,sl. Figure 8 illustrates the log-log regression 
for the datasets in the lower row of Figure 1. In addi- 
tion to the corresponding point estimate, we provide 
a 90% central interval estimate, using the paramet- 
ric bootstrap method as proposed by Davies and 
Hall (1999). 

Constantine and Hall (1994) argued that the bias 
of the variogram estimator increases with the cut- 
off L in the log-log regression, and Davies and Hall 
(1999) showed numerically that the mean squared 
error (MSE) of the estimator for a Gaussian pro- 
cess with powered exponential covariance is min- 
imized when L = 2. Zhu and Stein (2002) argued 
similarly in a spatial setting. These results are un- 
surprising and have intuitive support from the fact 
that the behavior of the theoretical variogram in an 
infinitesimal neighborhood of zero determines the 
fractal dimension of the Gaussian sample paths. We 
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Fig. 8. Log-log regression for the variogram estimator (13) and the datasets in the lower row of Figure 1. Only the two points 
marked with filled circles are used when fitting the regression line. 
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thus choose L = 2 in our implementation, resulting 
in the estimate 



(13) Dy;: 



llogV^(2/n 
2 



\ogV2{l/n) 



log 2 



As Kent and Wood (1997) suggested, generalized 
least squares rather than ordinary least squares re- 
gression could be employed, though the methods co- 
incide when L = 2. 

It is well known that the method of moments esti- 
mator (11) is nonrobust. It is therefore tempting to 
replace it by the highly robust variogram estimator 
proposed by Genton (1998), which is based on the 
robust estimator of scale of Rousseeuw and Croux 
(1993). We implemented an estimator of the fractal 
dimension that uses (12) with L = 2, but with the 
method of moments estimate (11) replaced by Gen- 
ton's highly robust variogram estimator. In a simu- 
lation setting, this estimator works well. However, it 
breaks down frequently when applied to real data, 
where it yields very limited, discrete sets of possible 
values for the estimate only. Upon further investi- 
gation, this stems from the ubiquitous discreteness 
of real-world data, under which the Rousseeuw and 
Croux (1993) estimator can fail, in ways just de- 
scribed. While discreteness is a general issue when 
estimating fractal dimension, the problem is exacer- 
bated by the use of this estimator. In this light, the 
next section investigates another approach to more 
robust estimators of fractal dimension. 

2.4 Variation Estimators 

We now discuss a generalization of the variogram 
estimator that is based on the variogram of order p 
of a stochastic process with stationary increments, 
namely, 

(14) -fp{t) = ^E\Xu-Xt+u\P. 

When p = 2, we recover the variogram (2), when p = 
1 the madogram, and when p=l/2 the rodogram 
(Bruno and Raspa, 1989; Emery, 2005; Bez and Ber- 
trand 2011). Standard arguments show that a Gaus- 
sian process with a variogram of the form (3) admits 
analogous expansions of the variogram of order p, in 
that 



as t — 7- 0, 

with fixed values of the fractal index, a E (0,2], 
/3 > 0, and a constant Cp> that satisfies 

\ 2/(ap) 



|(a+/3)p/2N 



IT 



C2- 



The fractal index, a, of the Gaussian process and 
the Hausdorff dimension, D, of its sample paths then 
admit the linear relationship (4). 

A natural generalization of the method of mo- 
ments variogram estimator (11) for time series or 
line transect data of the form (6) is the power vari- 
ation of order p, namely. 



(16) Vpil/n) 



1 



2{n-l) 



V IP 
^{i-l)/n\ ■ 



We then define the variation estimator of order p for 
the fractal dimension as 



-}S^(si-s)logVp{l/n) 



(17) 



. 1=1 



. 1=1 



where L > 2, s/ = log(//n) and s is the mean of 
Si,. . . ,sl. This definition nests the variogram, mado- 
gram and rodogram estimators, which arise when 
p = 2, 1 and 1/2, respectively. The general case, 
p > 0, has been studied by Coeurjolly (2001, 2008). 

For the same reasons as before, and supported by 
simulation experiments, we let L = 2 in our imple- 
mentation, so that 



(18) Dy., 



^ llogVp{2M-logVp{l/n) 



P 



log 2 



As an illustration, Figure 9 shows the log-log re- 
gression fit for the variation estimator of order p = 
1 and our example data. For instances of the use 
of the madogram estimator in the applied litera- 
ture, see Weissel, Pratson and Malinverno (1994) 
and Zaiser et al. (2004). 

A natural question then is for the choice of the 
power index p > 0. With the estimator depending on 
the relationship (4) between the fractal index, a, in 
the expansion (15) and the Hausdorff dimension, D, 
of the sample path, it is critically important to as- 
sess its validity when the assumption of Gaussian- 
ity is violated. In the standard case in which p = 2, 
Hall and Roy (1994) showed that, while the relation- 
ship (4) extends to some non-Gaussian processes, 
it fails easily. For example, it applies to marginally 
power-transformed Gaussian fields if and only if the 
transformation power exceeds 1/2. Other counterex- 
amples can be found in the work of Bruno and Raspa 
(1989) and Scheuerer (2010). 
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Fig. 9. Log-log regression for the madogram estimator [variation estimator (18) with power index p = IJ and the datasets in 
the lower row of Figure 1. Only the two points marked with filled circles are used when fitting the regression line. 



Bruno and Raspa (1989) and Bez and Bertrand 
(2011) applied the Lipschitz-Holder heuristics of 
Mandelbrot [(1977), page 304] to argue that the rela- 
tionship (4) is universal when p = \. While we agree 
that the relationship is particularly robust and in- 
clusive when p = 1, the Lipschitz-Holder heuristics, 
which connects the Lipschitz exponent of a sample 
path to its Hausdorff dimension, is tied to continuity. 
Thus, it can fail if the sample paths are sufficiently 
irregular. For instance, the sample paths of a binary 
stochastic process, which attains the values and 1 
only, have Hausdorff dimension D = 1. As the corre- 
sponding variogram (14) is independent of its order, 
we may refer to the common version as 7. If the ex- 
pected number of sample path jumps per unit time 
is finite, then 7 grows linearly at the coordinate ori- 
gin (Masry, 1972). In this case, the relationship (4) 
holds if, and only if, the common variogram, 7, is un- 
derstood to be of order p = 1. However, there are bi- 
nary processes whose variogram behaves like 7(t) = 
Odil"^) as t — )• 0, where < 7 < 1, and then the rela- 
tionship fails. Notwithstanding these examples, the 
argument of Bruno and Raspa (1989) and Bez and 
Bertrand (2011) is persuasive, and we maintain that 
the relationship (4) is particularly inclusive when 
p = 1. A natural conjecture is that if p = 1, then 
the relationship is valid if the process is ergodic (in 
a suitable sense) and the expected number of sam- 
ple path discontinuities per unit time is finite. Fur- 
thermore, it is worth noting that there are processes 
for which the madogram exists and the foregoing is 
satisfied, while second moments do not exist (Ehm, 
1981). 

The following interesting connection between the 
Hall-Wood estimator and the madogram estimator 
also suggests a special role of the power index p = l. 



For I a positive integer and j = 0, 1, — 1, let 

/ L(n-i)/«J 

A^'j\l/ n) = - ^ \X(^ii+j)/n - ^(i/+j-i)/n| • 
i=l 

Then A{1 / n) = A^^\l / n) and 

i>i(Vn) = i^i^l(^)(//n) 
i=o 

is, up to inessential constants, the mean of / dis- 
tinct copies of A{l/n). A comparison of the gen- 
eral forms (9) and (17), or the standard forms (10) 
and (18), of the Hall-Wood estimator with the mado- 
gram estimator suggests that the latter is a statisti- 
cally more efficient version of the Hall- Wood estima- 
tor. A similar, more tedious calculation can be used 
to argue heuristically that the box-count estimator 
has a bias, typically leading to lower estimates of 
the fractal dimension than the Hall-Wood and vari- 
ation estimators. For a confirmation in simulation 
studies, see Section 3.2. 

Here we restrict attention to a small initial study 
that assesses the efficiency and outlier resistance of 
variation estimators. Figures 10 and 11 show the 
root mean squared error (RMSE) of the variation 
estimator (18) from Gaussian sample paths in de- 
pendence on the power index, p. The curves are 
computed from 1,000 independent realizations with 
sample size n = 1,024, correspond to fixed values of 
the fractal index, a, and have their minima marked. 
Figure 10 concerns the ideal Gaussian model, where 
the estimator performs best for power indices of 
about p = 2, corresponding to the variogram esti- 
mator, similarly to the observations of Coeurjolly 
[(2001), page 1417]. Figure 11 shows that the RMSE 
can deteriorate considerably under a single additive 
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Fig. 10. Root mean squared error (RMSE) of the variation estimator (18) as a function of the power exponent, p, for 
Gaussian fractional Brownian motion (left) and Gaussian processes with powered exponential covanance (right). The scale 
parameter used is c — 1, and each RMSE is computed from 1,000 Monte Carlo replicates of Gaussian sample paths under the 
sampling scheme (6), where n = 1,024. The curves correspond to fixed values of the fractal index, a, and have their minima 
marked. 
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Fig. 11. Same as Figure 10, except that in each sample path a randomly placed observation is contaminated by an additive 
Gaussian outlier with standard deviation 0.1. 
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outlier, with the effect being stronger for smoother 
sample paths, that is, higher values of the fractal in- 
dex. The smaller the power index, the more outlier 
resistant the variation estimator. 

A possible variant of the variation estimator uses 
p- moments of higher increments, as proposed by Kent 
and Wood (1997) and Istas and Lang (1997) in the 
case p = 2. For example, one could turn to second 
differences, rather than first differences, and base 
a log-log regression on 



1 



n 



(19) 



2{n-2l) 

n—l 



V 

n I 5 



i=l 



rather than the standard variation (16). Also, if more 
than two points are used in the regression, the gen- 
eralized least squares method could be used in lieu of 
the ordinary least squares technique. However, there 
is no clear advantage in doing so in applied settings, 
in which the corresponding covariance structure is 
unknown and needs to be estimated as well. 

2.5 Spectral and Wavelet Estimators 

We now consider the semi-periodogram estima- 
tor of Chan, Hall and Poskitt (1995), which oper- 
ates in the frequency domain and is closely related 
to the spectral estimator of Dubuc et al. (1989a). 
The basis for this estimator is the well-known fact 
that the spectral density function for a stationary 
stochastic process with a second-order variogram of 
the form (3) decays like Iwl"""""^ as frequency |a;| — )• 
oo (Stein, 1999). For a stationary Gaussian process 
{Xt:t€ [0,1]}, Chan, Hall and Poskitt (1995) de- 
fined 

B{uj) = 2 [ Xtcos{uj[2t-l])dt 
Jo 

and called 

J{uj) = B{Lof 

the semi-periodogram. Under weak regularity con- 
ditions, the expected value of the semi-periodogram 
decays in the same way as the spectral density func- 
tion. Suppose now that we have = 2m + 1 obser- 
vations Xt at times t = i/{2m), where i = 0, 1, . . . , 2m. 
In this setting, Chan, Hall and Poskitt (1995) ap- 
proximated B{uj) by 



Biu;) 



1 



m 



V \ V 2m-l 

2 ^ 

i=l 



X, 



i/{2m) 



cos UJ 



m 



m 



and the semi-periodogram J{oo) by 
J{uj) = B{ujf. 

The semi-periodogram estimator of the fractal di- 
mension D is 



Dp 



(20) 




^{si - s)\ogJ{ui) 
1=1 

\ -1 




where w/ = 27r/, si = logw/ and s is the mean of si, 
. . . ,sl. The highest unaliased frequency (i.e., the Ny- 
quist frequency) is vrm, which is reflected by the fact 
that B['Km + 5) = B{'Km — 5) for any 5. This sug- 
gests setting L = [m/2j , but Chan, Hall and Poskitt 

(1995) recommended using L= [min{m/2, }J , 
which is less than \in/2\ for m > 34. As m grows, 
this rule thus has the effect of eliminating J{uj) at 
high frequencies in forming Dp , which at first seems 
counterintuitive, given that the underlying scaling 
law applies as uj increases. However, as uj approaches 
the Nyquist frequency, aliasing causes the expecta- 
tion of J{u}) to deviate markedly from the decay 
rate of t<;~"~^, leading to the need to eliminate high- 
frequency terms when forming Dp. Figure 12 shows 
an example of the log-log regression for the semi- 
periodogram estimator for datasets from Figure 1, 
where = 1,025 and hence L = 101 < m/2 = 256. 

The definition of B(uj) is similar in spirit to the 
so-called type-II discrete cosine transform (DCT-II); 
see Ahmed, Natarajan and Rao (1974) and Strang 
(1999) for background. Davies (2001) noted that this 
transform has some attractive properties when used 
as a basis for spectral analysis, so it is of interest 
to explore the DCT-II as a substitute for the semi- 
periodogram in estimating fractal dimension. Given 
time series data of the form (6), taking the definition 
of the DCT-II given by Gonzalez and Woods (2007) 
and adjusting it for our convention for the samples 
Xi/(2m) yields 

1/2 2m 



B{u) 



2m + 1 



i=0 



i/{2m) 



COS LO 



2i + l 
Am 



and J{uj) = B{lo)'^. Here the Nyquist frequency is 
27rm, as can be seen by noting that J(27rm + 6) = 
J(27rm — 6) for any 6. If we now let coi = 27rlm/(2m + 
1) with si = logw; and s being the mean of the 
s;'s as before, the log-log regression estimator (20) 
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Fig. 12. Log-log regression for the semi-periodogram estimator and the datasets in the lower row of Figure 1. Only the points 
marked with filled circles are used when fitting the regression line. 



can be applied with J(a;/) replaced by J{'jJi) and 

L= [min{2m,4ns }J . The DCT-II estimator uses 
approximately four times more points in the log-log 
regression than does the semi-periodogram estima- 
tor. For example, in Figure 13 we have L = 406, 
whereas we had L = 101 in Figure 12. 

The semi-periodogram estimator also serves as mo- 
tivation for a similar wavelet estimator, which is 
an adaptation of a weighted least squares estima- 
tor for the long memory parameter of a fractionally 
differenced process (Percival and Walden, 2000, Sec- 
tion 9.5). Given a time series of length n^, we com- 
pute its maximal overlap discrete wavelet transform 
(MODWT) out to level Jq = [log2(n^)J using reflec- 
tion boundary conditions; this can be done using the 
function modwt in the R package wavelets (Aldrich, 
2010). This MODWT yields Jq vectors of wavelet 
coefficients Wj, j = 1,..., Jq, each of which con- 
tains 2ns coefficients. The coefficients in the j'th vec- 
tor are associated with the scale Tj = 2^~^. The aver- 
age of these coefficients squared, that is, || Wj |p/2ns, 
provides an estimator of the wavelet variance v'^ (tj ) . 
This variance varies approximately as rj* for large Tj 
[Percival and Walden, 2000, equation (297b)], 



where a is the fractal index, from which the fractal 
dimension can be deduced; see Table 2. The scale tj 
corresponds to the band of frequencies ,7r/2^~^] 
The information that is captured by the semi- 
periodogram at high frequencies is thus captured in 
the wavelet variance at small scales tj. Since Chan, 
Hall and Poskitt (1995) eliminated certain high fre- 
quencies in their semi-periodogram estimator, this 
suggests using just the wavelet variances indexed 
by j = Jo, • ■ • , -^1, where Jq = max{l, [log2(ns)/3 - 
Ij}. Because the variance of the wavelet variance 
estimators depends upon tj, we replace the ordi- 
nary least squares estimator of the slope that is 
the basis for equation (20) with a weighted least 
squares estimator, say awL [Percival and Walden, 
2000, equation (376c)]. The corresponding estima- 
tor of the fractal dimension D is -DwL = 2 — |q;wl- 
Figure 14 gives an example of the wavelet estimator 
of D. Note that the estimator of D in the right-hand 
plot is -DwL = 2.05, which is greater than the upper 
limit D = 2 for the Hausdorff dimension of a curve. 
Sampling variability can cause this to happen on oc- 
casion with the other estimators also. In the simu- 
lation experiments reported below, the wavelet esti- 




FlG. 13. Log-log regression for the DCT-II estimator and the datasets in the lower row of Figure 1. Only the points marked 
with filled circles are used when fitting the regression line. 
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Fig. 14. Log-log regression for the wavelet estimator and the datasets in the lower row of Figure 1. Only the points marked 
with filled circles are used when fitting the regression line. 



mator proved to perform comparably to the DCT-II 
estimator, so we have chosen to drop the former and 
report only on the latter in what follows. 

3. PERFORMANCE ASSESSMENT: TIME 
SERIES AND LINE TRANSECT DATA 

We now turn to an evaluation of the various types 
of estimators, where we consider the large sample 
behavior under infill asymptotics and report on a fi- 
nite sample simulation study that assesses both ef- 
ficiency and robustness. 

3.1 Asymptotic Theory 

As noted, the fractal dimension refers to the prop- 
erties of a graph in a hypothetical limiting process 
that might exist if the scale of measurement were to 
become infinitely fine. Hence, estimators of fractal 
dimension are studied under infill asymptotics (Hall 
and Wood, 1993; Stein, 1999), in which the number 
of observations grows to infinity, whereas the un- 
derlying domain, namely the unit interval, remains 
fixed. We assume that time series or line transect 
data of the form (6) arise from a Gaussian process 
{Xt :t G [0, 1]} with a second-order structure of the 
type (3), where we let n grow without bounds. Typi- 
cally, the literature assumes stationarity, so that the 
process {Xt :t £ [0, 1]} has a covariance function of 
the form 

cj(t) = cj(0) - |ct|° + 0(|t|"+^) ast^O, 

where a £ (0,2), /? > and c > 0. The behavior of 
the bias, variance and mean squared error (MSE) 
of the estimators, and the corresponding types of 
limit distributions, then depend on the fractal in- 
dex a and on f3. Typically, the corresponding asymp- 
totic results carry over to Gaussian processes with 
stationary increments and a variogram or structure 
function of the form (3). 



For any Hall-Wood or variogram estimator D of 
the form (9) or (12) with a fixed value of the de- 
sign parameter L, the key results of Hall and Wood 
(1993) and Constantine and Hall (1994) are that 

MSE{D) 

(21) 

rO(n-i) + 0(n-2/3), ifO<a<f, 
= < 0(n-Mogn) + 0(n-2/5), ifa = |, 
io(n2"-4) + 0(n-2/3), if|<a<2, 

where in each case the first term corresponds to the 
variance, and the second term to the squared bias. 
If a < 3/2, then D has a normal limit; if a > 3/2, 
the limit is a Rosenblatt distribution as described by 
Taqqu (1975). In a recent far-reaching paper, Coeur- 
jolly (2008) showed that in the Gaussian case and for 
the variation estimator (17) with general power in- 
dex p> 0, the asymptotic behavior is still described 
by (21). Furthermore, the convergence rates are re- 
tained if the arithmetic mean in the definition of the 
power variation (16) is replaced by a trimmed mean, 
or by a convex combination of sample quantiles. 
While some of these results carry over to certain spe- 
cific non-Gaussian processes (Chan and Wood, 2004; 
Achard and Coeurjolly, 2010), the limiting distribu- 
tion theory is considerably richer then, and a general 
non-Gaussian theory remains lacking. 

It is interesting to observe the change in the asymp- 
totic rate of convergence at a = 3/2 for all these 
types of estimators. However, Kent and Wood (1997) 
showed that the variogram estimator achieves an 
MSE of order 

(22) MSE{D) = 0{n-^) + 0{n-^'^) 

for all a G (0, 2) if the design parameter satisfies 
L > 3 and the generalized least squares technique, 
rather than the ordinary least squares method, is 
used in the log-log regression fit, and/or second dif- 
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ferences of the form (19) are used. Similarly, Coem'- 
jolly (2008) demonstrated that the asymptotic rate 
of convergence for the variation estimator with gen- 
eral power index p > can be improved if second 
differences or related special types of increments are 
used. In finite sample simulation studies for the var- 
iogram estimator [p = 2), Kent and Wood (1997) 
did not find a clear-cut advantage in using general- 
ized least squares and/or second differences, and our 
own experience with variation estimators of diverse 
power indices is similar. As Kent and Wood (1997) 
argued, the likely cause is that, the closer a is to 2, 
the larger n must be before the asymptotic regime 
is reached. This behavior is in marked contrast to 
the case of spatial lattice data, to be discussed be- 
low. 

Chan, Hall and Poskitt (1995) developed asymp- 
totic theory for the semi-periodogram estimator, but 
the MSE decays at best at rate 0{n~^^^) in their 
results. The asymptotic scenario for the level cross- 
ing estimator in Feuerverger, Hall and Wood (1994) 
involves a bandwidth parameter and thus is not di- 
rectly comparable. 



3.2 Simulation Study: Gaussian Processes 

We now turn to a simulation study, in which we 
confirm and complement the foregoing asymptotic 
results in a Gaussian setting. In doing so, exact 
simulation is critically important (Chan and Wood, 
2000; Zhu and Stein (2002)), and we use the cir- 
culant embedding approach (Dietrich and Newsam 
1993; Wood and Chan, 1994; Stein, 2002; Gneit- 
ing et al., 2006) as implemented in the R package 
RandomFields (Schlather, 2001) to generate Gaus- 
sian sample paths, using the function GaussRF. The 
circulant embedding technique relies on the fast Fou- 
rier transform and is both exact and fast. 

Figure 15 shows log-log plots for the root mean 
squared error (RMSE) of the various types of es- 
timators in their dependence on the sample size n, 
computed from 1,000 independent trajectories of the 
form (6) from the corresponding stationary Gaus- 
sian process with a powered exponential covariance 
function, a{t) = exp(— |t|°^). The graphs are approx- 
imately linear, and their slopes show good agree- 
ment with the asymptotic laws in (21). Further- 
more, they confirm the aforementioned observation 
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Fig. 15. Root mean squared error (RMSE) of estimators of fractal dimension in dependence on the sample size, n, com- 
puted from Gaussian sample paths of the form (6) with powered exponential covariance function, a{t) = exp( — For each 
combination of a and n, the number of Monte Carlo replicates is 1,000. 
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Fig. 16. Boxplots for estimates of the fractal dimension from Gaussian sample paths of the form (6), where n = 1,024, 
with powered exponential covariance function, a{t) = exp( — and the fractal index, a, being equal to 0.4, 1-0 and 1.6, 
respectively. The corresponding true values of the fractal dimension, D, namely 1.8, 1.5 and 1.2, are shown as dashed lines. 
The number of Monte Carlo replicates is 500. 



of Kent and Wood (1997) that large values of the 
fractal index, a, require large sample sizes to reach 
the asymptotic regime. The variogram estimator gen- 
erally shows the lowest MSE, followed by the mado- 
gram and rodogram, and then the Hall- Wood, DCT- 
II, periodogram and box-count estimators. This rank- 
ing is retained under Gaussian processes with covari- 
ance functions from the Matern and Cauchy fami- 
lies, as well as for fractional Brownian motion, for 
all values of the fractal index a and all sufficiently 



large sample sizes, n. The use of variation estima- 
tors based on second differences, as defined in (19), 
typically does not yield lower RMSEs (results not 
shown) . 

Figures 16 and 17 show box- and scatterplots for 
the same types of estimators and the same class of 
Gaussian processes, where the sample size is n = 1,024 
Three groups of estimators can be distinguished, 
the first comprising the variogram and other two 
variation estimators along with the Hall- Wood esti- 
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Fig. 17. Scatterplot matrix for estimates of the fractal dimension from Gaussian sample paths of the form (6) with exponential 
covariance function, a{t) = exp( — and sample size n = 1,024. The true value of the fractal dimension is D — 1.5. The panels 
along the diagonal show histograms of the estimates, and those above the diagonal pairwise Pearson correlation coefficients. 
The number of Monte Carlo replicates is 500. 



mator, the second the spectral estimators, and the 
third the box-count estimator. The most efficient 
estimator is the variogram estimator, closely fol- 
lowed by the madogram estimator. As we have ar- 
gued theoretically before, the madogram estimator 
is a more efficient version of the Hall- Wood estima- 
tor, in that the estimators are strongly correlated, 
but the former is less dispersed. While the spectral 
estimators are less competitive, showing substan- 
tially higher dispersion than the variation estima- 
tors, the DCT-II estimator improves considerably 



on the periodogram estimator. The box-count es- 
timator generally shows a bias, with the estimates 
being too low. 

Figure 18 illustrates these results in a further ex- 
periment, in which we consider a Gaussian sam- 
ple path with the exponential covariance function, 
a{t) = exp(— and estimate the fractal dimension 
along sliding blocks of size 1,024. The corresponding 
estimates are plotted at the midpoint of the sliding 
block. It is clearly seen that the variogram and other 
variation estimators are the least dispersed, and that 
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Fig. 18. Estimates of the fractal dimension plotted at the midpoints of a sliding estimation window of size 1,025, for a Gaus- 
sian sample path of the form (6), where n — 10,000, with exponential covariance function, a{t) — exp( — The true fractal 
dimension, D = 1.5, is marked by the dashed line. The label on the horizontal axis, i, indicates the midpoint of the sliding 
estimation block, at i/10,000. 



the madogram estimator is a more efficient version 
of the Hall- Wood estimator. The spectral estimators 
are the most dispersed, with the DCT-II estimator 
outperforming the periodogram estimator, and the 
box-count estimator is biased. 

Thus far in this section, we have studied the ef- 
ficiency of the estimators under an ideal Gaussian 
process assumption. However, robustness against de- 
viations from Gaussianity is a critically important 
requirement on any practically useful estimator of 
the fractal dimension. In this light, we now expand 
our simulation study, and consider a situation in 
which Gaussian sample paths are contaminated by 
additive outliers. Specifically, given a sample path 
of the form (6), we let i be discrete uniform on 
{0, 1, . . . , n} and replace by -|- y, where y 
is normal with mean zero and standard deviation 
0.1 and independent of i. This process is repeated 
to obtain the desired number of outliers. 

Figure 19 shows RMSEs from such an experiment, 
using sample size n = 1,024, five additive outliers, 
and the powered exponential covariance function, 
a{t) = exp(— with values of the fractal index a 
that nearly span the full range from to 2. Not sur- 
prisingly, the results resemble those in Figures 10 
and 11, which considered variation estimators and 
the case of a single outlier only, and echo the find- 
ings of Achard and Coeurjolly (2010). Amongst the 
variation estimators considered here, the most out- 



lier resistant is the rodogram estimator {p = 1/2), 
but the box-count estimator, which performs poorly 
overall when there are no outliers, becomes compet- 
itive at the highest a considered (1.9). 

3.3 Discussion 

The foregoing results and arguments lead us to 
a recommendation for practitioners, in that we join 
Bruno and Raspa (1989) and Bez and Bertrand (2011) 
and call for the use of the madogram estimator, 
that is, the variation estimator with power index 
p= 1. The madogram estimator can be interpreted 
as a statistically superior version of the Hall- Wood 
estimator, is simultaneously more outlier resistant 
and more efficient than many of its competitors, 
and has strong intuitive appeal. Importantly, the 
critical relationship (4) between the fractal index 
of a stochastic process whose variogram of order p 
shows a behavior of the form (5) at the origin, and 
the fractal dimension of its sample paths, is appar- 
ently valid for a larger class of non-Gaussian pro- 
cesses when p = 1 than for certain other choices of p 
(in particular, p = 2), thereby justifying the use of 
the madogram estimator for both Gaussian and non- 
Gaussian stochastic processes. Its resistance to out- 
liers can be enhanced further if the arithmetic mean 
in the definition of the power variation (16) is re- 
placed by a trimmed mean, as proposed by Coeur- 
jolly (2008). 



20 



T. GNEITING, H. SEVCIKOVA, AND D. B. PERCIVAL 




Fig. 19. Root mean squared error (RMSE) of estimators of fractal dimension in dependence on the fractal index, a, for Gaus- 
sian sample paths of the form (6) with sample size n = 1,024 and powered exponential covariance function, a{t) = exp( — 
The panel on the left corresponds to the ideal Gaussian process setting; the panel on the right to a situation with five additive 
outliers in each sample path. The number of Monte Carlo replicates is 1,000. 



4. ESTIMATING THE FRACTAL DIMENSION 
OF SPATIAL DATA 

We now turn to estimators of the fractal dimen- 
sion of spatial data, as discussed by Dubuc et al. 
(1989b), Constantine and Hall (1994), Davies and 
Hall (1999), Chan and Wood (2000) and Zhu and 
Stein (2002), among other authors. Burrough (1981) 
noted a wealth of applications to landscape and other 
environmental data, with those of Rothrock and 
Thorndike (1980) on the underside of sea ice, and 
Goff and Jordan (1988) on the topography of the 
sea floor, being particularly interesting examples. 

From a probabilistic point of view, a natural ini- 
tial question is for the theoretical relationship be- 
tween the fractal dimension of a surface indexed 
in M^, and the fractal dimension of its sections or 
line transects. Assuming stationarity of the spatial 
random field and additional (weak) regularity condi- 
tions, Hall and Davies (1995) showed that the frac- 
tal dimensions along line transects are all identical 
to one another, except that in one special direction 
the dimension may be less than in all others. Very 
general results that do not depend on the stochas- 
tic process setting are available from the fundamen- 
tal work of Marstrand (1954). We join Davies and 
Hall (1999) in arguing that these results provide sub- 
stantial support for the use of fractal dimension as 



a canonical measure of surface roughness. In partic- 
ular, they allow us to estimate the fractal dimen- 
sion of a surface by adding 1 to any estimate of 
the fractal dimension of the corresponding line tran- 
sects. 

Technically, we focus discussion on the situation 
in which a spatial stochastic process, indexed by the 
unit square in M^, is sampled on a regular lattice, to 
yield a surface graph of the form 



i2 = 0,l,...,n^ C 



(23) 



Before reviewing estimators of fractal dimension stud- 
ied in the extant literature, and introducing new 
estimators, we propose a simple, unified notation. 
Specifically, for A: > we let 



S{k) = {{ii,i2,ji,j2) G {0,1,..., n} 



4 . 



k 



and denote the cardinality of this set by N{k). If 
N(k) > 1, we refer to A; as a relevant distance. The 
estimators in the subsequent Sections 4.1 to 4.2 then 
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take the form 



(24) 



5 = 2 - - 1 (sfe - s) fog Vp{k/n) I 



keJC 



where /C is a finite cohection of refovant distances, 
Sfc = fog(/c/n), s is the mean of {sk-k G /C}, and 
Vp{k/n) is a certain variation with general power 
index p> 0. 

Two-dimensional geometry allows for many op- 
tions in the choice of the distance representatives 
and the variation, and we restrict attention to the 
most plausible and best performing estimators in the 
literature, all of which are based on power variations 
(Davies and Hall, 1999; Chan and Wood, 2000; Zhu 
and Stein (2002)). As concerns the set fC of distance 
representatives, simulation experiments, experience 
in the line transect case, and the work of Chan and 
Wood (2000) and Zhu and Stein (2002) all suggest 
that a restriction to the smallest relevant distances 
only tends to lead to the best performance. In ad- 
dition to minimizing the bias of the estimator, this 
strategy keeps the computational complexity low as 
well. 

4.1 Isotropic Estimator 

Davies and Hall (1999) considered an estimator 
based on the isotropic empirical variogram, which 
we now generalize. For a relevant distance k, con- 
sider the variation 

VlSO;p(^/"') 

(25) 

■ \-^it/n,i2/n 



2N{k) 



S[k) 



with general power index p > 0. The isotropic es- 
timator -Diso;p with power index p then is defined 
by (24) with the set /C = {1, \/2, 2} of distance rep- 
resentatives and the variation V given by (25). Thus, 
we consider variations at horizontal and vertical 
spacings of one and two grid points (A; = 1 and k = 2), 
and a diagonal spacing of a single grid point (fc = 
V2), respectively. 

4.2 Filter Estimator 

Zhu and Stein (2002) studied a broad range of 
increment-based estimators, among which the "Fil- 
ter 1" estimator shows good performance. We gen- 
eralize by defining a filter estimator with general 
power index p> 0, rather than just p = 2 as in the 



work of Zhu and Stein (2002). Specifically, for a rel- 
evant distance > let 

^'^^^^^Z") = 2im 



(26) 



■ y ^j\-^ii/n,i2/n 
S(k) 



-2X 



(n+ii)/(2n),(i2+i2)/(2n) 
"T ^ji/n.j2/nl • 



The filter estimator D-p-p with power index p then 
is defined by (24) with the set K. = {2,2^2,4} of 
distance representatives and the variation V given 
by (26). This considers variations at horizontal and 
vertical spacings of one and two grid points {k = l 
and k = 2), and a diagonal spacing of a single grid 
point [k = V2), respectively. Hence, the filter esti- 
mator is t^he natural equivalent of the isotropic es- 
timator -Diso;p) but now using second differences, 
rather than first differences. 

4.3 Square Increment Estimator 

The square increment estimator is based on a pro- 
posal of Chan and Wood (2000), who restricted 
attention to quadratic variations. Here we define 
a square increment variation of general power index 
p> 0, namely. 



(27) 



2N{k) 

■ \-^ii/n,i2/n 
S{k) 



ii/n,j2/n 



Y J- Y IP 

^ji/n,i2/n ~r ^ji/n,j2/n\ 7 

where A; is a relevant distance. The square incre- 
ment estimator Dsi-p then is defined by (24) with 
the set fC = { ^2, 2^2} of dist ance representatives, 
corresponding to squares that have side widths of 
one and two grid points, and the variation V given 
by (27). 

4.4 Transect Estimators 

Finally, we consider two very simple estimators 
that are based on the variation estimator Dy-p with 
general power index p > of Section 2.4, or a vari- 
ant that uses second differences, as defined in equa- 
tion (19). In either case, a line transect estimate of 
the fractal dimension is computed for each row and 
each column in the grid. The transect- variation and 
transect-increment estimators -DtV;p and -DtI;p with 
power index p then add 1 to the median of the 2n 
corresponding line transect estimates. In the former 
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case, the line transect estimates are based on first 
differences, in the latter on second differences. 

5. PERFORMANCE ASSESSMENT: SPATIAL 
DATA 

We now assess the estimators by their large sam- 
ple behavior under infill asymptotic as well as in 
finite sample simulation studies, considering both 
efficiency and robustness. 

5.1 Asymptotic Theory 

Davies and Hall (1999), Chan and Wood (2000) and 
Zhu and Stein (2002) developed asymptotic theory 
for a very wide range of estimators of the form (24) 
that are based on variations with power index p = 2. 
Generally, their results apply under an infill asymp- 
totic scenario for sample paths of the form (23) from 
an intrinsically stationary Gaussian spatial process 
with fractal index a £ (0, 2) and a variogram that 
behaves like (3) at the origin. For estimators that 
are based on variations corresponding to a first dif- 
ference, the generic result is that 

MSE(5) 

(28) 

( 0(n-2) + 0{n-^^), if < a < 1, 

= < 0{n-^L{n)) + 0(n-4/3), if q = 1, 
[c>(n2"-4) + 0(n-4/3), ifl<a<2, 

where L is a function which is slowly varying at in- 
finity. If a < 1, then D has a normal limit, while, 
if Q > 1, the limit is related to a Rosenblatt distri- 
bution, with some of these results carrying over to 
certain specific non-Gaussian processes (Chan and 
Wood, 2004). However, if the variations correspond 
to a second difference, such as in the cases of the 
filter and square increment estimators, and/or the 
generalized least squares techniques, rather than the 
ordinary least squares method, is used, an improved 
asymptotic behavior, namely, 

(29) MSE{D) = 0{n~^) + Oin-^^) 

for all < a < 2, is achieved, with an associated 
limit distribution that is normal. For regularity con- 
ditions and further detail, we refer to the original 
work of Davies and Hall (1999), Chan and Wood 
(2000), Zhu and Stein (2002) and Chan and Wood 
(2004), which is impressive. While these authors re- 
stricted attention to quadratic variations with power 
index p = 2 only, we conjecture, based on the work of 
Guyon and Leon (1989), Istas and Lang (1997) and 
Barndorff-Nielsen, Corcuera and Podolskij (2009) 
on the power variations of Gaussian processes, that 



analogous results continue to hold under a general 
power index p> 0, similar to the line transect case 
studied by Coeurjolly (2008). 

As the amount of data in the sampling scheme (23) 
is about n^, the asymptotic rates of convergence 
in (28) and (29) conform with those in the time se- 
ries or line transect case, except that the transition 
from the classical to slower rates of convergence oc- 
curs already at a = 1 {ov D = 3/2), rather than at 
= 3/2 (or Z) = 7/4). Hence, these results suggest 
that there may be a higher benefit to using incre- 
ments that are based on second differences in the 
spatial case than in the time series or line transect 
case. 

5.2 Simulation Study: Gaussian Spatial 
Processes 

In the subsequent simulation study, we use state- 
of-the-art implementations of the circulant embed- 
ding method (Stein, 2002; Gneiting et al., 2006) to 
generate Gaussian sample surfaces. Figure 20 shows 
the root mean squared error (RMSE) of the estima- 
tors versus the (square root of the) sample size, n, 
using data of the form (23) from stationary spa- 
tial processes with powered exponential covariance 
function, a{t) =exp(— |t|"). The estimators use the 
traditional power index p = 2, and the number of 
Monte Carlo replicates is 1,000. 

The asymptotic rates of convergence in (28) 
and (29) suggest linear graphs with slope —1 on 
the logarithmic scale for the filter and square in- 
crement estimators and all values of the fractal in- 
dex, a. For the isotropic estimator, they suggest 
slope —1 for a < 1, slope —1/2 for a = 3/2, and 
slope —1/4 for a = 7/4. Our empirical results are 
in good agreement with the theoretical slopes, and 
attest to the superior performance of the filter es- 
timator in the ideal Gaussian process setting, as 
noted by Zhu and Stein (2002). Also, these and 
other simulation results lead us to conjecture that 
the transect-variation estimator behaves like (28), 
while the transect-increment estimator shares the 
favorable uniform asymptotic rate of convergence 
in (29). 

In Figure 21 we consider estimators with general 
power index p > 0, but fix n = 256 in the spatial 
sampling scheme (23). Again, we use the powered 
exponential covariance model, and the number of 
Monte Carlo replicates is 1,000, each comprising a to- 
tal of (n + 1)^ = 66,049 observations within the unit 
square. The left column shows the RMSE in the 
ideal Gaussian process setting, in which the filter. 
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Fig. 20. Root mean squared error (RMSE) for estimators of fractal dimension for data of the form (23) from spatial stochastic 
processes with powered exponential covariance, a{t) — exp(— versus the (square root of the) sample size, n. The power 
index used is p = 2, and for each combination of a and n, the number of Monte Carlo replicates is 1,000. 



square increment and transect-increment estimators 
perform best. Furthermore, the efficiency of these 
estimators depends only very little on the choice of 
the power index. 

Figure 21 also studies the behavior of the esti- 
mators in the presence of outliers. Specifically, the 
middle and right-hand columns show RMSEs in sit- 
uations in which the Gaussian sample paths have 
been contaminated by 10 and 20 additive outliers, 
respectively, in ways essentially identical to those 
described in Section 3.2. Note that the vertical scale 
differs from row to row, with the largest RMSEs cor- 
responding to the smoothest surfaces. The smoother 
the surface, that is, the larger the value of the fractal 
index, a, the more impact the outliers have on the 
estimators. Two of the estimators that dominate in 
the uncontaminated setting, namely the filter and 
the square increment estimators, are the least out- 
lier resistant, even though the outlier fraction is very 
small in our study, at 0.015 and 0.03 percent, respec- 
tively. The most resistant estimators are the transect 
estimators. 

These results allow for an interpretation in terms 
of breakdown points. While comprehensive formal 



definitions of breakdown points for dependent data 
have recently become available (Genton and Lucas, 
2003), it suffices here to take a heuristic point of 
view, and define the breakdown point of an estima- 
tor as the minimal fraction of contaminated data 
that can ruin an estimator. Evidently, the isotropic, 
filter and square increment estimators have break- 
down point zero. In contrast, the transect estimators 
have a positive breakdown point of about l/(2-^/rn) 
under the spatial sampling scheme (23), where m = n? 
is the approximate amount of data. To see this, note 
that each outlier affects the individual estimate for 
at most two of the 2n transects, from which the me- 
dian is formed. Thus, a transect estimator cannot 
be ruined, unless a fraction of at least (n/2)/n^, or 
l/(2-y/m), of the data are contaminated. 

5.3 Discussion 

While confirming extant theoretical and simula- 
tion results for spatial data, which suggest the use 
of variations that are based on second differences, 
the results of this section lead us to two novel in- 
sights. 
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Fig. 21. ijoot mean squared error of isotropic, filter, square increment, transect-variation and transect-increment estimators of fractal dimension in dependence 
on the power index, p. The underlying sample paths are of the form (23) with n — 256 from Gaussian spatial processes with powered exponential covariance, 
a(t) = exp( — \t\'^) . The columns correspond to situations with no outliers (left), 10 additive outliers (middle) and 20 additive outliers (right), respectively. The rows 
are for distinct values of the fractal index, namely a — 0.5, 1.0, 1.5 and 1.75. The number of Monte Carlo replicates is 1,000. 
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Thus far, the statistical hterature has restricted 
attention to variation-based estimators with power 
index p = 2. A first observation is that the efficiency 
of these estimators depends httle on the power in- 
dex, p > 0. Outher resistance and robustness argu- 
ments then suggest the use of smaller power indices, 
with Section 2.4 supporting the choice of p = 1. 

A second and potentially very surprising insight 
is the superior performance of the transect-variation 
and transect-increment estimators. These estimators 
have positive breakdown points and outperform the 
traditional estimators even under minimal devia- 
tions from the ideal Gaussian process setting. In- 
deed, in practice, we would expect much larger de- 
viations than in our simulation setting, in which the 
outliers had substantially lower variability than the 
process itself, and occurred at fractions of about 1 in 
6,000 and 1 in 3,000 only. As the transect-increment 
estimator appears to share the favorable uniform 
rate of convergence (29) with the best performing 
extant estimators, this suggests the availability of an 
estimator, namely the transect-increment estimator 
with power index p = l, that is both robust and ef- 
ficient. We believe that these are highly promising 
prospects that deserve further study. 

If spatial data are observed at irregularly scat- 
tered locations, rather than a regular grid, the only 
available estimator is a suitably modified version 
of the isotropic estimator (25), with S{k) now rep- 
resenting pairs of stations that are approximately 
a distance k apart. In such cases, the use of the 
power index p=l again seems prudent. 

6. DATA EXAMPLE: ARCTIC SEA-ICE 
PROFILES 

In this penultimate section, we estimate the frac- 
tal dimension of arctic sea-ice profiles based upon 
measurements of sea-ice draft (93% of thickness). 
These profiles can be regarded as line transects 
through the underwater surface of sea ice because 
they were collected using upward-looking sonars by 
submarines traveling under the sea ice with no ap- 
preciable deviations from a single direction and 
depth. The data used here were sampled one meter 
apart at a resolution of one meter by the U.S. Navy in 
August 1998 in the Arctic Ocean, and are available 
online from the National Snow and Ice Data Center 
at http : //nsidc . org/data/docs/noaa/g01360_ 
upward_looking_sonar. We examined six profiles of 
about 240 km total length (files sc98drft.002, 003 
and 005-008). These profiles are illustrated in Fig- 
ures 22 and 23 and show pronounced non-Gaussian 



features. Given the resolution, we join the extant 
literature in using fractal dimension to character- 
ize the surface roughness of the macroscopic, topo- 
graphic structure of sea ice, so that meters corre- 
spond to sufficiently small scales. 

We use a sliding estimation window of width n = 
1,024 meters, move these blocks along the profiles 
in increments of 10 meters, and estimate fractal di- 
mension for each of them, using variation estima- 
tors with power indices p = 2 (variogram), 1 (mado- 
gram) and 1/2 (rodogram) along with the Hall- 
Wood estimator. Thus, for each method there are 
23,303 blocks in total, with Figure 23 showing exam- 
ples of profiles and dimension estimates. Figures 24 
and 25 show the corresponding boxplots, histograms 
and pairwise scatter diagrams, composited over all 
blocks. While all four methods result in similar esti- 
mates, the strongest correlation is, not surprisingly, 
between the Hall-Wood and madogram estimators, 
with the latter being our estimator of choice. 

Overall, the Arctic sea-ice profiles have Hausdorff 
dimension of about 1.3 along line transects, and thus 
of about 2.3 spatially. These estimates are at the 
lower end of results reported in the literature. For 
instance. Bishop and Chellis (1989) argued that pro- 
files of ice keels have fractal dimensions ranging from 
1.2 to 1.7, while Connors, Levine and Shell (1990) 
reported estimates of about 1.4 and 1.6 for first year 
and multiyear ice segments, respectively. Further re- 
sults and background information can be found in 
the works of Rothrock and Thorndike (1980) and 
Goff et al. (1995). 

7. CONCLUDING REMARKS 

In closing this review, we return to the opening 
quote of Lenny Smith [(2007), page 115] that pre- 
cedes the abstract of our paper. Indeed, healthy skep- 
ticism about dimension estimates based on real- world 
data is well justified, in that estimates of fractal di- 
mension depend on the availability of data at suit- 
able scales and face issues of discreteness, and the ef- 
fects of measurement error might also need to be dis- 
entangled. Notwithstanding these issues, estimates 
of fractal dimension can serve as informative de- 
scriptors of surface roughness. 

In the case of time series or line transect data, 
we recommend the use of the madogram estimator, 
that is, the variation estimator with power index 
p = l. In the case of spatial lattice data, transect es- 
timators based on madogram estimators along rows 
and columns show high promise. These recommen- 
dations echo observations of Bruno and Raspa (1989) 




and Bez and Bertrand (2011), who argued that the 
critical relationship (4) between the fractal index of 
a stochastic process whose variogram of order p> 
shows a behavior of the form (5) at the origin, and 
the fractal dimension of its sample paths, is partic- 
ularly robust and inclusive when p= 1. We encour- 
age work toward rigorous results in these directions, 
including both variograms and their equivalents for 
general types of line transect and spatial increments. 

Furthermore, we call for the development of large 
sample theory for variation estimators of general 
power index, including but not limited to our pre- 
ferred choice of p = 1. While for time series this 
has been achieved in the far-reaching recent work 
of Coeurjolly (2008), the case of spatial data re- 
mains open. As noted in Section 5.1, we believe 
that many of the results in the extant default case 
p = 2 can be carried over to a general power in- 
dex, p > 0. In doing so, the results of Guyon and 
Leon (1989), Istas and Lang (1997) and Barndorff- 
Nielsen, Corcuera and Podolskij (2009) on the power 
variations of Gaussian processes provide tools that 
can be applied in concert with the methodology put 
forth in an impressive strand of asymptotic litera- 



ture for p = 2, which includes the work of Hall and 
Wood (1993), Constantine and Hah (1994), Kent 
and Wood (1997), Davies and Hah (1999), Chan 
and Wood (2000) and Zhu and Stein (2002), among 
others. Of particular interest is our conjecture in the 
spatial setting of Section 5.1, according to which the 
transect-variation approach allows for estimators of 
the fractal dimension that are simultaneously highly 
efficient and robust. This approach can be paired 
with the use of trimmed means or linear combina- 
tions of sample quantiles in defining power varia- 
tions, as suggested by Coeurjolly (2008) for time 
series, or with the use of bipower and multipower 
variations (Barndorff-Nielsen et al. 2006; Barndorff- 
Nielsen, Corcuera and Podolskij, 2011). While the 
case of Gaussian processes appears to be challenging 
yet tractable, a general asymptotic theory for non- 
Gaussian processes remains elusive, despite the re- 
cent progress by Chan and Wood (2004) and Achard 
and Coeurjolly (2010). 

To our knowledge, Bayesian methods for estimat- 
ing fractal dimension have not been explored yet, ex- 
cept that the method of Handcock and Stein (1993) 
could be applied in a parametric context. Physical 
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Fig. 23. Selected blocks of width 1,024 meters from the sea-ice profiles in Figure 
of fractal dimension. 




'flOlOO 



mhowk i.n 



ISOM IDZnO 25400 ZEm SHOD 

along with the corresponding estimates 



insight can drive the choice of the prior, and the de- 
velopment of Bayesian madogram estimators might 
be a promising option. 

Multifractional Brownian motion (Peltier and Levy 
Vehel, 1995; Benassi, Jaffard and Roux, 1997; Herbin, 



2006) is a generalization of the classical fractional 
Brownian motion, where the fractal dimension is al- 
lowed to vary along the sample paths. Similarly, the 
nonstationary Gaussian random fields described by 
Anderes and Stein (2011) allow for location-depend- 
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Fig. 24. Boxplots for variation estimates with power indices p = 2 (variogram), 1 (madogram) and 1/2 (rodogram) along 
with Hall-Wood estimates of fractal dimension for the blocks m the ice profile data. See the text for details. 
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Fig. 25. Scatterplot matrix for variation estimates with power indices p = 2 (variogram), 1 (madogram) and 1/2 (rodogram) 
along with Hall-Wood estimates of fractal dimension for the blocks m the ice profile data. See the text for details. 



ent, local Hausdorff dimensions, which need to be 
estimated as functions, rather than a single num- 
ber. In this setting, the estimators considered in our 
paper can be used as building blocks for the more 
complex estimators needed to handle these nonsta- 
tionary processes, as studied by Benassi, Cohen and 
Istas (1998), Ayache and Levy Vehel (2004) and 
Coeurjolly (2005), among others. For an applied per- 
spective, see Gagnon, Lovejoy and Schertzer (2006). 

To close on a practical note, we have developed an 
R package, called f ractaldim, that implements the 
estimators of fractal dimension discussed in this pa- 
per (Sevcikova, Gneiting and Percival, 2011). It has 
the ability to compute estimates for a single dataset, 
or a series of estimates along sliding blocks, as in 
our data example. Multiple estimates can be conve- 
niently bundled into a single call, and the default 
arguments correspond to the recommendations in 
this paper. 

For example, to generate four log-log plots of the 
type shown in Figure 5, a few lines of code suffice: 

par (mf row=c (2,2)) 

series <- GaussRF(x=c (0, 1 , 1/1024) , 
model= ' stable ' , gridtriple=TRUE, 
param=c(mean=0, variance=l, 
scale=l, kappa=l)) 



methods <- cChallwood' , 

'periodogram' , 'variogram' , 

'madogram' ) ) 
D <- f d. estimate(series , 

method=methods , plot . loglog=T, 

plot . allpoints=T) 

Here, the function GaussRF from the RandomFields 
package is used to generate a Gaussian sample path 
with exponential covariance and fractal index a = 
1. The desired estimation methods are passed to 
the estimation function, fd. estimate, as charac- 
ter strings as in our example. Alternatively, each 
method can be given as a list with a name element 
determining the method, and any additional ele- 
ments specifying parameters. For example, an entry 
list (name=' variation' , p.index=l) is equivalent 
to 'madogram'. 

To obtain interval estimates with a parametric 
bootstrap method, as proposed by Davies and Hall 
(1999) and shown in many of our figures, the follow- 
ing call suffices, using the simulated series from 
above, boot . it bootstrap replicates, and the mado- 
gram estimator: 

D <- fd.estimate(series, 

methods= ' madogram ' ) 
D.boot <- rep(NA, boot.it) 
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alpha <- 4 - 2*D$fd 
for (i in l:boot.it) 
{ 

boot. series <- GaussRF( 

x=c(0, 1, 1/1024) , model= ' stable ' , 

gridtriple=TRUE , 

param=c (mean=0 , variance=l, 

scale=D$scale~ (-1/alpha) , 

kappa=alpha) ) 
D.boot[i] <- fd. estimate ( 

boot . series , 

methods= ' madogram ' ) $f d 

> 

Confidence intervals can be obtained from the ar- 
ray D.boot in the usual way. Additional functions 
for variation estimators in the line transect case are 
available within the dvfBm package (Coeurjolly, 
2009). 
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